Event Driven Langevin simulations of Hard Spheres 
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The blossoming of interest in colloids and nano-particles has given renewed impulse to the study 
of hard-body systems. In particular, hard spheres have become a real test system for theories and 
experiments. It is therefore necessary to study the complex dynamics of such systems in presence 
of a solvent; disregarding hydrodynamic interactions, the simplest model is the Langevin equation. 
Unfortunately, standard algorithms for the numerical integration of the Langevin equation require 
that interactions are slowly varying during an integration time-step. This in not the case for hard- 
body systems, where there is no clearcut between the correlation time of the noise and the time-scale 
of the interactions. Starting first from a splitting of the Fokker-PIank operator associated with the 
Langevin dynamics, and then from an approximation of the two-body Green's function, we introduce 
and test two new algorithms for the simulation of the Langevin dynamics of hard-spheres. 

PACS numbers: 05.10.-a,05.40. Jc,05.10.Gg 



I. INTRODUCTION 

Hard spheres (HS) are a reference system for structural 
and dynamical theories of fluids [11,0, but idealized: the 
infinitely steep potential is essentially a way of capturing 
the effects of steric interactions. On the atomic or the 
molecular scale two body interactions are mostly mod- 
elled by Lennard- Jones or Coulumb potentials; experi- 
ments on colloids shift the length scales of interest up to 
roughly 1 nm to 1000 nm where objects can behave as 
hard bodies and are still small enough to exhibit ther- 
mal or Brownian motion in a solvent. Dynamical light 
scattering [s', has already provided a rich collection of 
data for such systems, encouraging a considerable effort 
in understanding the dynamics; the possibility of follow- 
ing single particle trajectories via confocal microscopy 
of latex particle [5] has allowed a direct view on an ex- 
perimental realization of HS systems and their dynamics 

The simplest model of a suspension of neutral parti- 
cles is to consider a system of HS in an ideal solvent with 
no hydrodynamic interactions; real suspensions are of- 
ten described in terms of their deviations from such ideal 
system. This is the most interesting model for theoreti- 
cians and many results have been derived: the two body 
case (and hence the low density case) has been solved ex- 
actly [Si, 9] , while at moderate and high packing fractions 
various Enskog-like [Hi or Mode Coupling theories 
[l^ [Tsl have been applied to understand the dynamics. 
While hydrodynamic interactions (HI) are well under- 
stood at low particle densities, much less is known at 
high densities, and theories often proceed by disregarding 
them fl3] . As an example, theories regarding glass transi- 
tion often disregard HI effect, like in the case of the Mode 
Coupling theory for Brownian hard spheres [Tsi il6j or 
Brownian hard discs in shear fiow[l7|. 

Non-HI simulations therefore have their place in test- 
ing such theories, and in circumventing the huge effort 
needed to simulate HI [ll-[23. 



In order to validate non-HI theories for HSs it is nec- 
essary to use computer simulations, as only a qualita- 
tive agreement is to be expected among non-HI theories 
and data for real suspension. Standard simulation meth- 
ods for Brownian dynamic like the well-known Ermak- 
McCammon [2^ require continuous potentials; to cir- 
cumvent such problem several algorithms have been in- 
troduced with various degrees of justification 23-26] for 
the over-damped dynamics; only recently it has been rec- 
ognized that in the case of hard interactions such simu- 
lations are better performed by event-driven (ED) codes 
[27I - I29I . We introduce two new ED algorithms that go 
beyond the over-damped approximation and allow for the 
simulation of the full Brownian dynamics of HSs. 



II. IVIETHODS 

We consider a system of N HSs governed by the 
Langevin equation 



dtVi = 



(1) 



for the positions and the velocities v^; here 7 is the 
friction constant, = —ra"^ dyJJ the acceleration, m is 
the mass of the HSs, U is the potential energy and m^j 
are the zero-mean random forces due to the solvent. We 
assume that such random forces are delta correlated and 
satisfy the fluctuation-dissipation theorem 
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In the case of continuous interactions, it is possible to 
define stochastic Taylor expansions [30l | ; correspondingly, 
integration schemes of the fc-th order with errors of order 
{/Si)^ in the time-step Ai can be introduced 3l|. In the 
case of hard-body interactions, all the standard machin- 
ery of stochastic calculus breaks down due to the singular 
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nature of the interaction potential and new methods must 
be developed. 

We consider the Fokker-Plank equation associated to 
the SDE HI) (Kramers' equation [3a |) 

dtW = LKW (3) 

where W (r, v, t) is the probability distribution function 
(PDF) for the positions r — {rj} and the velocities 
V — {vi} of the particles, — kBT/m relates to the 
temperature and 

La' = 7 (^v • V + vl,dl) - (v + a • av) (4) 

is the Kramer operator. Integrating the SDE ([IJ for 
a finite time-step At corresponds to extracting a con- 
figuration |i-*+A*j v*"*"^*} according to the probability 
gL^f Atj _ V - V*). 

III. SPLITTED BROWNIAN DYNAMICS 

To obtain a numerical approximation, a powerful ap- 
proach is to split the evolution operator e^'^^^ in a prod- 
uct e^^^* « JJgQiLiAt exactly-integrable operators L.^ 

, , i 

[33| ensuring that that the decomposition is positive (i.e. 
all fli > 0) [34]. Therefore, to each splitting corresponds 
an algorithm in which in a single time-step At the oper- 
ators e"*'"''^* are applied in sequence. 

We first choose to split Lik into the reversible 
(or streaming) operator hrev — ~ {v ■ dr + a - d^) 
and the irreversible (or collision) operator L^^r = 
7 (9v • V -I- vff^d^) [35]; we indicate the corresponding al- 
gorithm as Splitted Brownian Dynamics (SBD). 

The operator Jjrev is the Liouvillian associated to the 
Hamiltonian H = mv ■ v/2 + U. In the case of step 
potentials, the associated reversible equation of motion 
can be integrated via event-driven molecular dynamics 
(EDMD) 36] with a precision limited only by the nu- 
merical round-off errors; therefore the propagator e^rcvi^t 
can be implemented with extreme accuracy. 

The operator L^^r corresponds to the interaction with 
the bath; the associated SDE is dt'v = —7V + ^ can be 
exactly integrated giving an explicit formula for the evo- 
lution v*+^* = e^-'-^'v*: 

<T = e-^^'vU + v'^,\(l-e-27A*)r (5) 

where F is a unitary Gaussian random variable and a G 
{x,y,z}. 

The algorithm for the single SBD time-step 
gLrc^AtgLirrAt consists therefore in an EDMD sim- 
ulation [s^ of length At followed by a thermalization 
of the velocities according to eq.®. We notice that the 
error is at most quadratic (as can be checked via Taylor 
expansion (J^^^vAt^Ur^At ^ gL^c At _^ q ^^^2-^ 

regards only in the dynamics; in fact, SBD is equivalent 
(upon identifying the angle a mixing reversible and 



irreversible evolution with cos (a) = e-T(*-*')) to the 
Generalized Hybrid Monte Carlo |37i and therefore ex- 
plores the canonical ensemble as long as the propagation 
steps g^r-c^At ^^ijirrAt ^g^^ exactly implemented (as in 
our case). 

It is therefore of interest to give some physical bounds 
on the magnitude of the feasible time-step At. First, we 
notice that for At — 00 the dynamics reduces to and 
MD simulations where velocities are extracted each At 
from a Maxwellian; therefore if the time-step is much 
bigger than the average inter-particle collision time, re- 
sults of classical MD are to be expected. Accordingly, 
we find that for big At the algorithm overestimates the 
diffusion coefficient (fig. [T]); this is to be expected as 
the mean free path (in absence of collisions) of a parti- 
cle is of order VthAt instead of ^~^VthVAi. Second, the 
magnitude of At is naturally bounded the damping time 
T = ; therefore the SBD is not well indicated for 
simulations in the over-damped limit j/m — >■ 00. Ac- 
cordingly, we find that SBD overestimates diffusion co- 
efficients for At > 7~^ (figlj); it is therefore necessary 
to develop an alternative approach for the simulation of 
systems with high damping. 

IV. APPROXIMATE GREEN'S FUNCTION 
DYNAMICS 

It has been shown in that the over-damped limit 
of eq.([T]) can be simulated efficiently using ED codes (2ft]. 
The algorithm relies on considering time steps At small 
enough so that mostly binary collisions are relevant, i.e. 
the average displacement should be less than the aver- 
age inter-particle separation. Moreover, average displace- 
ment should be smaller than the HSs' radii in order to 
map the interaction of two nearby HSs in the problem 
of a random walk near a reflective wall. Under such ap- 
proximations, the true two-body stochastic dynamics for 
over-damped Brownian HSs can be implemented by algo- 
rithm of 29] in which each step consists in predicting the 
displacements Ax of the HSs via the free propagator, in- 
troducing fictive velocities v-'^ — Ax/ At, and performing 
an EDMD with such fictive velocities during t and t + At. 
We extend such approach to the general Brownian case. 

First, we need to predict the positions of the HSs after 
a time-step At according to their free propagation, i.e. 
the solution of eq.([T]) with no interaction (a = 0): 

r v*+^* = V* -f- A^ + AvR , , 

I j.t+At ^ r-t + Ar + Atr ^ ' 

The particle displacements contain both systematic parts 
Av = (e >"* — 1) V*, Ar — 7~"'^(1 — e"'*'*) v* and stochas- 
tic displacements. The stochastic displacements Av^^, 
Ar^ are zero-mean correlated gaussian variables with 
variances (Av|.) = m-^kBT{l - e^^T*)^ (Ar|j) = 
7-im-ifci3r [2i-7-i (3 + 4e-T*-|-e-27t)] and cross- 
correlation (Arj^Avfl) = -f-'^m-^kBT (1 - e'^^f [H]. 



3 



If we consider a time-step such, that the average dis- 
placement is less than the average inter-particle separa- 
tion, we can consider only the corrections due to two- 
body interactions. In the limit of small At, a couple of 
HSs will interact only when they start from nearby posi- 
tions. In particular, if j^^vthVAt a, i.e. the average 
free displacement is much smaller than the diameter a 
of the HSs, the dynamics of two particles A and B can 
be approximated as the Langevin dynamics of a point 
particle at a distance (r^ — r^) (1 — a/ ||r^ — rsH) from 
a flat wall. It is possible to solve such problem with a 
straightforward generalization of the image method ap- 
plied in 29]. In fact, the solution given by the free 
particle Green's function plus an image particle with 
a reflected velocity beyond the reflective wall (fig. [5]) 
correctly satisfies the zero-current boundary condition 
^'ilwaii ~ 0, where n is the normal to the wall and 
j(r, i) = J vW{r,v,t)dv is the probability current for 
the position. 

Such solution can be implemented exactly by pre- 
dicting the new positions and velocities jV*"*"^* 
according to eq.®, defining Active velocities = 
^j.t+At _ j.t^ I ^-f- performing an EDMD simulation 
with such Active velocities during At; if a collision hap- 
pens, the component of the relative velocity normal to 
the contact point must be reflected for both the flctive 
and the predicted velocities v*"*""^'. We indicate such 
algorithm as the approximate Green's function dynamics 
(AGD). In the over-damped limit, the prediction of the 
velocities and positions decorrelates and the algorithm 
correctly reduces to the over-damped case of [1^ . 

As for the SBD algorithm, it can be proven that the 
AGD scheme respects detailed balance and ergodicity 
and therefore explores the correct ensemble for HSs; 
hence, errors are again only in dynamic quantities. At 
difference with SBD, we have no analytic estimate for 
the error; nevertheless, we expect that the the mean- 
free path in absence of collisions 'y~^vthVAt must be 
smaller than the radius of the HSs in order to satisfy the 
flat-wall approximation, and must be smaller than the 
average inter-particle distance in order to avoid multiple 
collisions (hence higher than two-body effects) during At. 

In order to check that the behaviour of AGS is driven 
just by geometrical considerations, we have simulated HS 
systems at different 7 and varying the time-step At in 
the range [l0~2, 10"] (reduced units). At difference with 
SBD where diffusion can vary even by a order of mag- 
nitude in such a At range, the values of D measured 
from AGD vary a few percent over the range and long 
simulations are been necessary to have enough statistics 
to detect the behaviour of D that would otherwise look 
flat. In fig. Owe show that the measured diffusion coeffi- 
cient D versus the AGS simulation time-step displays a 
plateau (i.e. fluctuations become much smaller than 1%) 
already for At < 0.1 regardless of 7 and 0. 



V. CONCLUSIONS 

Hard spheres, and in general hard body systems in 
suspension, have become a realistic model due to the 
developments of experimental techniques for the inves- 
tigation of colloidal systems and nano-particles; yet the 
dynamics of such systems is hard to simulate via the 
standard Brownian dynamics algorithms. In fact, classi- 
cal continuous-time algorithms fail due to instantaneous 
character of the interactions; we have shown instead how 
it is possible to simulate the full Langevin dynamics of 
Hard Spheres. 

First, we have shown how the simplest splitting of the 
stochastic evolution operator (a technique often referred 
to as " Trotterization" from Trotter's seminal work[39|) 
allows to write an algorithm (the Splitted Brownian Dy- 
namics - SBD). The SBD algorithm becomes inefficient of 
high viscosities but via the operator-splitting technique 
could easily take account for the interaction with exter- 
nal fields or with the presence of fiuxes (like shear) in the 
surrounding fiuid. 

Second, we have shown how by considering the two 
body dynamics of Brownian Hard Spheres it is possi- 
ble do develop an algorithm (the Approximate Green's 
function Dynamics AGD) that overcomes such problem 
and works equally well for a wide range of packing frac- 
tions and viscosities. To develop the AGD algorithm, 
we have solved the problem of the Langevin dynamics 
dfV = —^v + ^ oi a point particle in presence of a reflec- 
tive wall by extending the classical Image Method solu- 
tion for the over-damped Brownian dynamics dtx = 77 of 
a point particle in presence of a reflective wall (here ^, 
T] are noises). The AGD algorithm is Event Driven and 
considers fictive collisions between Hard Spheres. While 
it should possible to take into account the polydispersity 
of a system by considering also effective masses in the 
fictive collisions as in [2^, including shear or external 
fields in the AGD algorithm looks more complicated as 
it would require the solution of the particle - reflective 
wall problem with external fields/shear. 

Both SBD and AGD simulations explore the canon- 
ical ensemble for Hard Spheres and therefore repro- 
duce the correct equilibrium thermodynamics. They be- 
long to the class of Asynchronous Event-Driven Particle 
Algorithms |40|] and can be easily implemented by adapt- 
ing exis ting codes for ED dynamics [36j or Brownian Dy- 
namics \4V\ of Hard Spheres. 
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FIG. 1. Effect of the damping coefficient 7 on tlie size of 
tlie simulation step At (all quantities in reduced units). The 
diffusion coefficient D from simulations is plotted versus the 
time-step size At for various 7's. As expected, the system 
approaches the MD value for diffusion regardless of 7 for 
At 00. The "true" value of D is obtained for At — > . 
We observe at small At's a plateau in the D vs At plot for 
At < 7~^, signalling that the "true" value of D is approached. 
Results are presented for packing fraction — 0.30; a com- 
pletely analogous behaviour is found at a low packing fraction 
(l> ~ 0-10 and an high packing fraction (j> = 0.45. 
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FIG. 2. A two body problem for hard spheres can be mapped 
into the problem of a point particle interacting with a larger 
sphere. When particles are very near, the problem further 
simplifies to the interaction of a Langevin particle with a re- 
flective flat wall, whose solution can be derived by applying 
the Image Method to the Langevin equation. In fact, the 
Green function must zero inside the wall and must satisfy 
the no-flux boundary conditions at the wall. Combining the 
free Green function of the particle in its initial position and 
the free Green function of its image (with the normal-to-the- 
wall component of the initial velocity reflected) satisfles both 
Kramer' equations and reflective boundary conditions giving 
the correct solution. 



FIG. 3. Effects of packing fraction (j) (upper panel) and of 
the damping coefficient 7 (lower panel) on the time-step At 
for the AGF algorithm. All quantities in reduced units; thick 
lines are just a guide for the eye. Diffusion D is calculated 
averaging over 10 independent trajectories for 2000 particle 
systems; simulations are long at least 10 times the structural 
correlation time. In the upper panel, results are shown for cj> — 
0.10,0.30,0.45 at fixed damping 7 — 10. In the lower panel, 
results are shown for 7 = 1, 10, 100 at fixed packing fraction 
(j> = 0.30. Notice that the estimated diffusion coefficient D 
has a small relative variation in the wide range of dampings 
7s and packing fractions (j)s analysed. As a rule of thumb, to 
estimate D with an accuracy much smaller than 1% time-step 
of order At ~ 0.1 are already enough. 



